Read data into R
Several functions read data from The Climate Analyzer into R. Below are examples of how to use these functions to import annual, monthly, and daily temperature and precipitation data.
GHCN Co-op weather stations
First, we need to get a list of GHCN Co-op stations hosted on The
Climate Analyzer. The stations() function retrieves these
data. In this example we are going to use the type = "COOP"
argument to get a list of Co-op stations.
ghcn_stations <- stations(type = "COOP")
dplyr::glimpse(ghcn_stations)
#> Rows: 656
#> Columns: 9
#> $ station_name <chr> "Aberdeen Experimental Station", "Alta 1NNW", "Antelope…
#> $ id <chr> "USC00100010", "480140", "350197", "470349", "100470", …
#> $ type <chr> "COOP", "COOP", "COOP", "COOP", "COOP", "COOP", "COOP",…
#> $ lat <dbl> 42.95300, 43.77270, 44.81972, 46.56660, 44.04240, 32.33…
#> $ lon <dbl> -112.8250, -111.0340, -120.7533, -90.9660, -111.2740, -…
#> $ elev_m <dbl> 1342, 1962, 924, 198, 1589, 1151, 305, 200, 200, 1263, …
#> $ years_avail <chr> "1914;1915;1916;1917;1918;1919;1920;1921;1922;1923;1924…
#> $ years_segments <chr> "1914 - 2025", "1909 - 1940;1948 - 2025", "1924 - 1943;…
#> $ data_url <chr> "ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/all/USC001…There are 656 GHCN Co-op station hosted on The Climate Analyzer. Lets filter this down to a more focused area, like Canyonlands National Park (NP), Utah.
For this example we are going to use the arcgislayers
package to read in the National Park Service (NPS) park boundaries from
ArcGIS Online. The URL below is a publicly available REST endpoint
hosted by the National Park Service. These data need to be transformed
to EPSG:4326 (WGS 84) so they align with the GHCN Co-op station
coordinates. Next, we are going to use the sf package to
clip the GHCN Co-op stations to Canyonlands NP. Finally, we will include
the Hans Flat Ranger Station Co-op station just outside the park
boundary.
# URL for NPS boundaries on ArcGIS Online
rest_url = "https://services1.arcgis.com/fBc8EJBxQRMcHlei/arcgis/rest/services/NPS_Land_Resources_Division_Boundary_and_Tract_Data_Service/FeatureServer/2"
# Read feature layer and filter to Canyonlands NP
cany <- arcgislayers::arc_read(rest_url) |>
janitor::clean_names() |>
sf::st_make_valid() |>
dplyr::filter(unit_code == "CANY") |>
sf::st_transform(crs = "EPSG:4326") # Transform data to WGS 84
# Convert to sf object
ghcn_sf <- sf::st_as_sf(ghcn_stations, coords = c("lon", "lat"),
crs = "EPSG:4326") |>
sf::st_make_valid()
# Clip GHCN Co-op stations in Canyonlands NP
cany_ghcn <- sf::st_intersection(ghcn_sf, cany) |>
dplyr::select(-tidyselect::any_of(colnames(cany))) |>
# Include the Hans Flat RS station just outside the park boundary
dplyr::bind_rows(
dplyr::filter(ghcn_sf, station_name == "Hans Flat RS")
)# View data
dplyr::glimpse(cany_ghcn)
#> Rows: 5
#> Columns: 8
#> $ station_name <chr> "Canyonlands The Neck", "Canyonlands The Needle", "Gray…
#> $ id <chr> "421163", "421168", "XXX3", "XXX4", "423600"
#> $ type <chr> "COOP", "COOP", "COOP", "COOP", "COOP"
#> $ elev_m <dbl> 1790, 1495, 0, 0, 2012
#> $ years_avail <chr> "1965;1966;1967;1968;1969;1970;1971;1972;1973;1974;1975…
#> $ years_segments <chr> "1965 - 2025", "1965 - 2025", "1980 - 2025", "1980 - 20…
#> $ data_url <chr> "ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/all/USC004…
#> $ geometry <POINT [°]> POINT (-109.821 38.46), POINT (-109.759 38.167), POINT …Now let’s make a quick map using the maptiles and
tmap packages to visualize the GHCN Co-op stations in
Canyonlands National Park.
# Get base tiles
aoa_tiles <- maptiles::get_tiles(cany, provider = "Esri.WorldTerrain",
zoom = 12)
#> |---------|---------|---------|---------|=========================================
# Map
tmap::tm_shape(aoa_tiles, unit = "mi") + tmap::tm_rgb() +
tmap::tm_shape(cany) + tmap::tm_polygons(fill = 'darkseagreen2',
col = 'darkgreen',
fill_alpha = 0.7) +
tmap::tm_shape(cany_ghcn) + tmap::tm_symbols(shape = 20, fill = 'blue') +
tmap::tm_compass(type = "arrow", position = c("left", "top"),
text.size = 0.5, bg.color = 'white', bg.alpha = 0.5) +
tmap::tm_scalebar(breaks = c(0, 5, 10, 15, 20), position = c(0.3, 0.05),
text.size = 0.4, bg.color = 'white', bg.alpha = 0.5) +
tmap::tm_add_legend(type = "polygons", labels = "Canyonlands NP",
fill = "darkseagreen2", col = 'darkgreen',
title = "Legend") +
tmap::tm_add_legend(type = "symbols", labels = "Co-op Stations",
fill = "blue", size = 0.7) +
tmap::tm_layout(main.title.size = 1, legend.position = c(0.01, 0.15),
legend.text.size = 0.5, legend.title.size = 0.7,
legend.frame = TRUE, legend.bg.color = 'white',
frame = TRUE) +
tmap::tm_title("Canyonlands NP GHCN Co-op Stations")
Finding your station_id
Unfortunately, there isn’t a good way to read the
station_id into R. Sorry about that! :/
To get the station_id, go to The Climate Analyzer and select the
station you want to download data for and get the station id from the
URL. For this vignette we are going to pull data from the GHCN Co-op
station at the Island in the Sky Visitor Center in Canyonlands National
Park . This station is named “Canyonlands The Neck” on The Climate
Analyzer. The station_id is the last section of the URL:
canyonlands_theneck.
Importing normalized temperature and precipipation data
The normals() function imports NOAA calculated
temperature and precipitation normals for a specified reference period.
In this example, we are going to use the 1981-2010 reference period.
Other options available are 1971-2000 and 1991-2020.
isky_normals <- normals(station_id = "canyonlands_theneck",
ref_period = "1981-2010")
dplyr::glimpse(isky_normals)
#> Rows: 39
#> Columns: 5
#> $ station_id <chr> "canyonlands_theneck", "canyonlands_theneck", "canyonlands_…
#> $ id <chr> "421163", "421163", "421163", "421163", "421163", "421163",…
#> $ element <chr> "PRCP", "PRCP", "PRCP", "PRCP", "PRCP", "PRCP", "PRCP", "PR…
#> $ month <fct> Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec,…
#> $ value <dbl> 0.64, 0.55, 0.79, 0.71, 0.70, 0.44, 0.89, 1.11, 1.09, 1.30,…Import annual temperature and precipitation data
The import_annual() function imports annual temperature
and precipitation data for a specified station and time period. In this
example we are going to import all the data (1965 to 2025) for the
Island in the Sky Co-op station.
isky_annual <- import_annual("canyonlands_theneck", 1965, 2025)
dplyr::glimpse(isky_annual)
#> Rows: 62
#> Columns: 7
#> $ year <dbl> 1965, 1966, 1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 19…
#> $ prcp_cy <dbl> NA, 8.76, 8.66, 6.03, 12.38, 7.22, 7.45, 9.61, NA, NA, 10.07, …
#> $ prcp_wy <dbl> NA, 10.66, 8.54, 7.08, 11.96, 7.37, 5.22, 6.39, 12.94, NA, 11.…
#> $ tmax_cy <dbl> NA, 64.98, NA, 62.36, 63.73, 63.69, 62.94, 64.67, NA, NA, 62.9…
#> $ tmax_wy <dbl> NA, NA, 63.91, NA, 64.04, 63.25, 63.77, 64.99, 60.60, NA, 61.7…
#> $ tmin_cy <dbl> NA, 41.87, 41.97, 40.61, 42.55, 42.25, 41.60, 43.29, NA, NA, 4…
#> $ tmin_wy <dbl> NA, 42.50, 41.94, 41.15, 42.47, 41.88, 42.28, 43.28, 39.34, NA…Now let’s use the annual_figure() function to plot
annual maximum temperature (TMAX), minimum temperature (TMIN), and
precipitation (PRCP) data along with the 1981-2010 normals.
# Maximum Temperature (TMAX) ----
# Get annual maximum temperature normals
ann_tmax_normals <- dplyr::filter(isky_normals,
element == "TMAX" & month == "Annual")
# Plot annual maximum temperature data
my_title <- "Annual Maximum Temperature (1965-2025)"
annual_figure(x_var = isky_annual$year, y_var = isky_annual$tmax_cy,
normal = ann_tmax_normals$value, reference_period = "1981-2010",
area_color = "coral", line_color = "coral4",
my_title = my_title, my_ylab = "Temperature (°F)")
# Minimum Temperature (TMIN) ----
# Get annual minimum temperature normals
ann_tmin_normals <- dplyr::filter(isky_normals,
element == "TMIN" & month == "Annual")
# Plot annual minimum temperature data
my_title <- "Annual Minimum Temperature (1965-2025)"
annual_figure(x_var = isky_annual$year, y_var = isky_annual$tmin_cy,
normal = ann_tmin_normals$value, reference_period = "1981-2010",
area_color = "lightskyblue1", line_color = "dodgerblue2",
my_title = my_title, my_ylab = "Temperature (°F)")
# Precipitation (PRCP) ----
# Get annual precipitation normals
ann_prcp_normals <- dplyr::filter(isky_normals,
element == "PRCP" & month == "Annual")
# Plot annual precipitation data
my_title <- "Annual Precipitation (1965-2025)"
annual_figure(x_var = isky_annual$year, y_var = isky_annual$prcp_cy,
normal = ann_prcp_normals$value, reference_period = "1981-2010",
area_color = "palegreen1", line_color = "forestgreen",
my_title = my_title, my_ylab = "Precipitation (inches)")
Import monthly temperature and precipitation data
Next, we are going to use the import_monthly() function
to import monthly temperature and precipitation data for the Island in
the Sky Co-op station from 1965 to 2025. After importing the data, we
are going to create a water year (Oct-Sep) variable and a factor
variable for month names so we can plot the data by water year using the
append_water_yr_mth() function.
isky_monthly <- import_monthly("canyonlands_theneck", 1965, 2025) |>
append_water_yr_mth()
dplyr::glimpse(isky_monthly)
#> Rows: 732
#> Columns: 10
#> $ year <dbl> 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965…
#> $ month <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7…
#> $ prcp <dbl> NA, NA, NA, NA, NA, NA, 2.21, 0.68, 0.70, NA, NA, 0.91, 0.…
#> $ tmax <dbl> NA, NA, NA, NA, NA, NA, 90.61, 89.00, 74.18, NA, NA, 39.00…
#> $ tmin <dbl> NA, NA, NA, NA, NA, NA, 63.36, 61.04, 49.54, NA, NA, 23.84…
#> $ water_month <int> 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1…
#> $ month_name <fct> Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec…
#> $ season <chr> "Wi", "Wi", "Sp", "Sp", "Sp", "Su", "Su", "Su", "Fa", "Fa"…
#> $ water_yr <dbl> 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1966…
#> $ water_mth <fct> Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec…We can plot monthly maximum temperature (TMAX), minimum temperature
(TMIN), and precipitation (PRCP) data using the
monthly_figure() function. We are going to use the
cowplot package to plot the three figures in a single
column in one figure.
plot_title <- cowplot::ggdraw() +
cowplot::draw_label("Water Year Monthly Climate Data (1965-2025)",
fontface = 'bold', size = 16, x = 0.5, y = 0.5)
# Monthly TMAX figure
mth_tmax <- monthly_figure(
x_var = isky_monthly$water_mth, y_var = isky_monthly$tmax,
year_group = isky_monthly$water_yr, my_year = 2025, line_color = "coral4",
my_title = "A. Monthly Maximum Temperature",
my_ylab = "Temperature (°F)"
)
# Monthly TMIN figure
mth_tmin <- monthly_figure(
x_var = isky_monthly$water_mth, y_var = isky_monthly$tmin,
year_group = isky_monthly$water_yr, my_year = 2025, line_color = "dodgerblue2",
my_title = "B. Monthly Minimum Temperature",
my_ylab = "Temperature (°F)"
)
# Monthly PRCP figure
mth_prcp <- monthly_figure(
x_var = isky_monthly$water_mth, y_var = isky_monthly$prcp,
year_group = isky_monthly$water_yr, my_year = 2025, line_color = "forestgreen",
my_title = "C. Monthly Precipitation",
my_ylab = "Precipitation (inches)"
)
# Combine figures
cowplot::plot_grid(plot_title, mth_tmax, mth_tmin, mth_prcp,
labels = c(), ncol = 1, align = "v", vjust = 1, scale = 1,
rel_heights = c(0.1, 1, 1, 1))
Import daily temperature and precipitation data
Next, we are going to use the import_daily() function to
import daily temperature and precipitation data for the Island in the
Sky Co-op station from 1964 to 2025. Again, we are going to to add the
water year (Oct-Sep) and water month variables using the
append_water_yr_mth() function.
isky_daily <- import_daily("canyonlands_theneck", 1964, 2025) |>
dplyr::mutate(date = as.Date(date, format = "%m/%d/%Y"),
year = lubridate::year(date),
month = lubridate::month(date)) |>
append_water_yr_mth()
dplyr::glimpse(isky_daily)
#> Rows: 22,280
#> Columns: 12
#> $ date <date> 1965-01-01, 1965-01-02, 1965-01-03, 1965-01-04, 1965-01…
#> $ prcp_in <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
#> $ tmax_f <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
#> $ tmin_f <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
#> $ snow_depth_in <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
#> $ year <dbl> 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 19…
#> $ month <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ water_month <int> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,…
#> $ month_name <fct> Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, J…
#> $ season <chr> "Wi", "Wi", "Wi", "Wi", "Wi", "Wi", "Wi", "Wi", "Wi", "W…
#> $ water_yr <dbl> 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 19…
#> $ water_mth <fct> Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, J…I haven’t used the daily data much yet, so I don’t have any plotting
functions for daily data. However, you can use ggplot2 to
create your own plots. Here’s a quick example of how to plot daily
maximum temperature (TMAX) data for water year 2025.
# Filter data for water year 2025
isky_daily_2025 <- dplyr::filter(isky_daily, water_yr == 2025)
# Plot daily maximum temperature data
ggplot2::ggplot(data = isky_daily_2025, aes(x = date, y = tmax_f)) +
ggplot2::geom_line(color = "coral4") +
ggplot2::labs(title = "Daily Maximum Temperature, Water Year 2025",
x = "Date", y = "Temperature (°F)") +
ggplot2::theme_minimal()
Import departure data
The import_departure() function imports monthly
temperature and precipitation departure data for a specified station and
time period. In this example we are going to import monthly departure
data for the Island in the Sky Co-op station from 1965 to 2025 and add
the water year (Oct-Sep) and water month variables using the
append_water_yr_mth() function..
isky_depart <- import_departure("canyonlands_theneck", 1965, 2025) |>
append_water_yr_mth()
dplyr::glimpse(isky_depart)
#> Rows: 732
#> Columns: 10
#> $ year <dbl> 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965…
#> $ month <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7…
#> $ prcp_pctavg <dbl> NA, NA, NA, NA, NA, NA, 228.10, 87.94, 80.08, NA, NA, 166.…
#> $ tmax_depart <dbl> NA, NA, NA, NA, NA, NA, 0.11, 1.00, -4.32, NA, NA, 1.20, -…
#> $ tmin_depart <dbl> NA, NA, NA, NA, NA, NA, -1.74, -2.26, -5.16, NA, NA, 2.84,…
#> $ water_month <int> 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1…
#> $ month_name <fct> Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec…
#> $ season <chr> "Wi", "Wi", "Sp", "Sp", "Sp", "Su", "Su", "Su", "Fa", "Fa"…
#> $ water_yr <dbl> 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1966…
#> $ water_mth <fct> Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec…Now we can plot monthly temperature (TMAX & TMIN) and
precipitation (PRCP) departure data using the
departure_figure() function.
plot_title <- cowplot::ggdraw() +
cowplot::draw_label("Water Year Monthly Climate Departures (1965-2025)",
fontface = 'bold', size = 16, x = 0.5, y = 0.5)
# Departure TMAX figure
depart_tmax <- departure_figure(
x_var = isky_depart$water_mth, y_var = isky_depart$tmax_depart,
year_group = isky_depart$water_yr, y_intercept = 0, my_year = 2025,
line_color = "coral4",
my_title = "A. Monthly Maximum Temperature Departure",
my_ylab = "Departure (°F)"
)
# Departure TMIN figure
depart_tmin <- departure_figure(
x_var = isky_depart$water_mth, y_var = isky_depart$tmin_depart,
year_group = isky_depart$water_yr, y_intercept = 0, my_year = 2025,
line_color = "dodgerblue2",
my_title = "B. Monthly Minimum Temperature Departure",
my_ylab = "Departure (°F)"
)
# Departure PRCP figure
depart_prcp <- departure_figure(
x_var = isky_depart$water_mth, y_var = isky_depart$prcp_pctavg,
year_group = isky_depart$water_yr, y_intercept = 100, my_year = 2025,
line_color = "forestgreen",
my_title = "C. Monthly Precipitation Departure",
my_ylab = "Departure (%)"
)
# Combine figures
cowplot::plot_grid(plot_title, depart_tmax, depart_tmin, depart_prcp,
labels = c(), ncol = 1, align = "v", vjust = 1, scale = 1,
rel_heights = c(0.1, 1, 1, 1))